//
// This file is part of cpp-lib, a library of C++ functions and classes.
// Cpp-lib is copyright (C) 2004 and onwards, The CPL Interest Group.
// 
// This library is free software; you can redistribute it and/or modify it
// under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation; either version 2.1 of the License, or (at
// your option) any later version.
// 
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
// General Public License for more details.
// 
// You should have received a copy of the GNU Lesser General Public License
// along with this library (see the file COPYING); if not, write to the
// Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
// 02111-1307  USA
//
// $Id$
//


#include <iostream>
#include <sstream>
#include <numeric>

#include <cmath>

#include "cpp-lib/matrix.h"
#include "cpp-lib/matrix-io.h"
#include "cpp-lib/banded_matrix.h"
#include "cpp-lib/math-util.h"
#include "cpp-lib/container-util.h"
#include "cpp-lib/util.h"

using namespace cpl::math;

typedef matrix< double , fixed< 3 , 1 > > vec_3 ;
typedef matrix< double , fixed< 3 , 3 > > mat_3 ;


double init[] =
  {
     1 , -1 , 0 ,
     1 ,  1 , 0 ,
     0 ,  0 , 1
  } ;


mat_3 const O = initialize( init ) ;

int main() {

  debug_output( std::cout , lower_bandwidth( O ) ) ;
  debug_output( std::cout , upper_bandwidth( O ) ) ;
  debug_output( std::cout , upper_bandwidth( O + O ) ) ;

  // Doesn't work because the preprocessor doesn't know about C++
  // template syntax!!!
  // debug_output( std::cout , identity< double , 3 , 3 >() ) ;

  matrix< double , banded > B( 10 , 10 , 2 , 3 ) ;
  B = 4711. * identity< double , 10 , 10 >() ;
  
  debug_output( std::cout , lower_bandwidth( B ) ) ;
  debug_output( std::cout , upper_bandwidth( B ) ) ;
  debug_output( std::cout , lower_bandwidth( transpose( B ) ) ) ;
  debug_output( std::cout , upper_bandwidth( transpose( B ) ) ) ;

  debug_output( std::cout , B ) ;
  
  debug_output( std::cout , lower_bandwidth( B + B ) ) ;
  debug_output( std::cout , upper_bandwidth( B + B ) ) ;
  debug_output( std::cout , lower_bandwidth( B * B ) ) ;
  debug_output( std::cout , upper_bandwidth( B * B ) ) ;

  std::cout << O ;
  std::cout << transpose( O ) ;
  std::cout << O * transpose( O ) ;
  std::cout << O * identity< double , 3 , 3 >() ;


  std::cout << identity< double , 10 , 10 >() ;

  mat_3 OO = 2. * O ;
  std::cout << OO ;

  vec_3 x = column_vector( 2. ,  2. , 4. ) , 
        y = column_vector( 0. , -4. , 2. ) ;

  debug_output( std::cout , x ) ;
  debug_output( std::cout , y ) ;
  debug_output( std::cout , x + y ) ;

  debug_output( std::cout , inner_product( x , y ) ) ;

  debug_output( std::cout , transpose( y ) ) ;

  debug_output( std::cout , upper_bandwidth( x ) ) ;
  debug_output( std::cout , lower_bandwidth( x ) ) ;
  debug_output( std::cout , upper_bandwidth( transpose( y ) ) ) ;
  debug_output( std::cout , lower_bandwidth( transpose( y ) ) ) ;

  debug_output( std::cout , x * transpose( y ) ) ;

  matrix< double , banded > M( 3 , 3 , 1 , 1 ) ;
  double const phi = 10 * pi / 180 ;
  M( 1 , 1 ) = std::cos( phi ) ; M( 1 , 2 ) = -std::sin( phi ) ;
  M( 2 , 1 ) = std::sin( phi ) ; M( 2 , 2 ) =  std::cos( phi ) ;
  M( 3 , 3 ) = 1 ;

  debug_output( std::cout , M     ) ;
  debug_output( std::cout , M * M ) ;
  
  std::cout << "should be 2: " ;
  debug_output( std::cout , lower_bandwidth( M * M ) ) ;
  debug_output( std::cout , upper_bandwidth( M * M ) ) ;

  {

    std::ostringstream os ;

    os.precision( 16 ) ;
    os << M ;
    matrix< double , banded > MM( 3 , 3 , 1 , 1 ) ;
    std::istringstream is( os.str() ) ;
    is >> MM ;

    std::cout << "should be small (at most about " 
              << std::numeric_limits< double >::epsilon()
	      << "): "
	      << max_abs_norm( M - MM ) 
	      << '\n' ;

  }


  std::cout << "sequence M^k x (M rotates by 10 degrees about z axis)\n" ;
  while( 1 ) { y = M * x ; x = y ; debug_output( std::cout , x ) ; }

}
